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1. Introduction 

Global fluctuations of space-averaged order parameters in many statistical systems 
possess generally a non-gaussian behavior in their critical and non-mean field regime. 
One well-known example is given by the two-dimensional XY- model (2D-XY), whose 
magnetic fluctuations (specifically the absolute value of the global spin vector) can be 
studied within the spin- wave approximation on the low-temperature critical line [TJ [2j [3] , 
far enough from the Kosterlitz-Thouless transition [H El E] • One interesting feature and 
universal character of the resulting rescaled probability distribution function (PDF) is 
that it does not depend on the critical exponents, temperature and system size [7j. 
Moreover the exponential and double exponential asymptotic fall-off behaviors for 
negative and positive deviations respectively, up to some corrective terms, identify this 
PDF to a general Gumbel function parametrized by a non-integer. Gumbel functions 
describe the extrema distribution of a set of independent and identically distributed 
Gaussian variables. They are parametrized by an integer which is equal to one for the 
lowest value, two for the second lowest values and so on. The connection with the 
2D-XY-model is based on the existence of low excitations modes, the Goldstone modes, 
which tend to destroy the quasi-long range order at low temperature. In Fourier space, 
the modes of value q, are decoupled and have a mass linearly proportional to |q|, and 
therefore the low soft modes are the predominant contributions to the PDF, giving 
rise to a strongly non-Gaussian behavior and Gumbel limiting distribution with non- 
integer parameter tt/2. Simple one dimensional and self-critical models such as interface 
distribution or 1/ f noise [8] give also a Gumbel distribution with parameter equal to one 
and therefore strongly non-Gaussian. Gumbel distributions with generalized parameter 
are also common in one-dimensional non-equilibrium and correlated systems where the 
distribution can be cast into a product of single probabilities with independent but non- 
identically distributed variables [HI [H)] . In all these example, the Fourier modes have 
mass proportional to |q|, therefore leading to the same limiting distribution dominated 
by the soft modes. This function is constructed by isolating the non-singular part of the 
fluctuations when the system size is taken to infinity, usually this is done by rescaling the 
order parameter using the mean value and variance only. Self-organized systems [HI [12] 
and network applications are also important systems in the way that they often lead 
to exact mathematical results [131 EH]; f° r example the distribution of connectivity 
in random graphs [15J, width distribution of interfaces [33], or size distributions in 
replicative phenomena [T71 [18]. Generally, one expects to approach the critical point 
from the low temperature region where the distribution for the order parameter m 
(for example the magnetization per site in a spin system) can be generally written as 
P(m, L) = L^l v Q{mL^l v 1 L, hL 135 ^), where £ is the correlation length, h an external 
field conjugated to the order parameter m, L the typical system size and Q the limiting 
or asymptotic PDF independent of L. The critical exponents (/3, 5, u) depend on the 
specific model. For the classical 2D-XY model, where m is defined as the modulus of 
the total spin vector divided by the number of sites L 2 , this asymptotic distribution can 
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be very closely written as [3] 




Q(9)^e W (-^W 92/2e - a °J + 8n,/^9) (1) 

witha = l/24 + 7/(47r)-ln(47r)/(47r)-lnnr=i 1 -exp(-27r£;) /(2tt) ~ 0.113514, 
and constant g-i ~ 3.867 x 10~ 3 . Here 9 = (m— < m >)/o~, with a 2 the variance, 
is the standard rescaled parameter for defining correctly the non-singular limit. The 
asymptotic regimes when \9\ is large are important in defining the class of function 
where the PDF belongs to. For independent and identically distributed fluctuations, 
the PDF falls into the Gaussian class. Otherwise, the PDF usually falls into class of 
functions characteristic of generalized extreme value statistics, for example Gumbel, 
Weibull or Frechet [19J. They are however few correlated systems where the asymptotic 
PDF of global order parameter can be evaluated analytically [201 E] arid where the 
effect of the correlations can be measured. In this paper we focus essentially on the 
critical regime of a simple correlated spin-model, the well-known spherical model in 
dimension 2 < d < 4 [21J, where fluctuations of the order parameter (magnetization 
and its absolute value) can be studied exactly in the critical regime. This will give some 
physical insight on universal features belonging to critical correlated systems. 



2. Definition of the Model 

We consider the spherical model [21] defined by a set of N = L d scalar spins 
— oo < Si < +oo coupled together with nearest-neighbor ferromagnetic coupling 
on a d- dimensional lattice at temperature T = 1/(3. The spins are constraint by the 
condition £V Sf = N, and the Hamiltonian is given by 



Hp — —77 /] JijSiSj, (2) 

with = JJ2a^i,rj±e a where e Q=lj ... d are the unit spacing vectors of the lattice. 
For example, in d = 2, ei = x and e.2 = y- Following the standard techniques, 
see [211 E2] , the partition function is written using a Dirac delta function representation 
to impose the constraint on the spins 



Z N (h) 



ds 
2^ 



] J dSi exp 



-(3H + hJ2S l + {is + a) In- S i 



(3) 



where a is any arbitrary real scalar which does not change the value of the integral 
over s. The successive derivatives of Z^(h) with respect to h give the moments of the 
magnetization m := ^YlnSi/N which we want to study. We will also consider another 
partition function by replacing the field term h Si by h Si\, and from which we 
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Figure 1. Surface plot of the function $>(z,h) (vertical axis) in the complex plane, 
showing the location of the unique saddle point. Here h = K = 1 and N = 10 3 in 
d = 3. 



can derive the successive moments of the absolute value of the magnetization \m\. In 
general we can consider both fields and study the different moments from Z^(h, h). The 
integration over the Si can be done using Fourier transform Si = ■^'^2qSqe lcl ' Ti , with 
S-q = Sq, q = 27r(/ci, • • • ,kd)/L, and fcj = 0, • • • , L — 1. In this case the interaction term 
can be diagonalized directly such as 

- E J ij S i S j = E Sq^-q E C0S 

i,j q a 

with K := j3J. After Fourier transformation, the constraint contribution becomes 
simply 



(is + a) 



E % = E E SqSq,e^^ = (is + a) E S q S_ 



(4) 



The other linear terms are equal to h ^ Si = hyN So, and h \ Si\ = hy iV j *So | - 
We can finally express the partition function as an integral over a product of decoupled 
Fourier modes 



Z N (h,h) 



ds 
2^ 



+oo 



dS Yl dSqdS-q exp 
q^O 



— K E SqS-q \z + d — E COS Qa 
q V a / 
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hVNSo + hVN \S \ + NK(z + d) 



(5) 



where K(z + d) = is + a. a is chosen to avoid the integral divergence for the zero 
mode. The prime symbol on the product means that only half of the Fourier modes 



are considered, using the symmetry q — > — q. Decomposing Sq 



it/q in real 



and imaginary parts, so that dSqdS-q = dx q dy n , S'qS'-q = x q + y^, we obtain after 
integration over Fourier modes the integral expression of the partition function 



Z N (h, h) 



ds ( ix \ N / 2 

2^ Km) exp 



- ln ( z + d - Yl cos ^ ) + N 

q \ a / 



h 2 + h 2 

AKz 



+ NK(z + d) 



X 



y/2 



2 cosh + e -^/(«*)erf ( VN ( h Jl h) ) + e^/(2^) erf ( ^ N ( h + h ) 



2Kz 



2VKz 



2VKz 



From this expression, we can define the intensive free energy h, h) for the 
system such that 



/da f IT \ N/2 
2^\2k) >/2ejq>[JV$(*,M) 

with, for each of the two order parameters considered in this paper 



(6) 



$(z, h,0) 

Hz Ah) 



h 2 



q \ a 



cos g a 



4Kz 



(7) 



1 + erf 



Nh 



2VKz 



It is clear that $ is not an even function of h for the absolute value of magnetization. 
As TV is large, we can study the saddle points of $ in the complex plane, and obtain the 
main contributions to the free energy and partition function. 



3. Saddle point analysis and scaling relations in the critical region 

In Figure HJ we represented $ for some arbitrary values of parameters, which shows the 
location of the unique saddle point as the minimum of $ along the real z-axis, which is 
also a maximum along the real s-axis or imaginary z-axis. In particular the equation of 
the local saddle point satisfies for example 



dHz,hA 1 v 1 1 ^ + K = ( 8 ) 

dz 2N ^ z + d-Ycos q a 2Nz AKz 2 ' 1 ' 

q^O ^ a HOi 

The last equation gives the implicit saddle point solution z(h) as function of field 
h. The critical point K = K c is defined by K c := lim7v-s.oo jfi X] q ^o(^ — S Q cos la) 1 ='■ 
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|<7i, after taking first the limit iV — > oo then h = and finally z — 0. This 
quantity is well defined for d > 2. In two dimensions for example, this quantity 
is pertinent for characterizing the XY model and diverges with the system size as 
gi ~ (1/271") log(CiV) [3] with the exact value 



C = exp|| + 21og (— ) +27-41ogJJ[l-exp(-27m)]| ~ 1.8456 

n=l 

The evaluation of this constant requires a precise analysis of the discrete sum over 
the Fourier modes which is specific to the geometry of the lattice. By extension it is 
useful to define the dimension-dependent quantities 



q^O V ^ a HOL * q^O 

which scale with L in the large size limit, and where G q is the Green function. 
Indeed, in the interval 2 < d < 4, we can replace the sum over modes q by an integral 
and show that 



1 rQmax T d 

g '^Lwr G ^ ldq ' (10) 

where q max is some cut-off dependent of the lattice step. The integrand is dominant 
when q is small, precisely G q oc g~ 2 , then g p m L d<yl ~ p ^> [q d ~ 2p ]l™ t a *L- ^ P = 1' ^ ne hitegral 
is convergent for d > 2, and g\ is finite, as well as the critical temperature. Otherwise, 
for p > 2, d — 2p is strictly negative if d < 4, and the integration part [q d ~ 2p ] 9 c ™ t a / x L is 
dominated by L 2p ~ d which is divergent. Therefore g p scales like £ J d ( 1 -p)+ 2 P- d = L-p( d - 2 ) 
and goes to zero like 



9 P >2 = ~g p L-^ + 0(L- d ^), ~g p := Km ^ G q , (11) 

q^O 

where g p are finite and positive constants. The corrective terms 0(L~ d ( p ~^) are 
always negligible since 2p > d, which holds when 2 < d < 4. The case d = 2 is 
particular since all g p are finite, except for g\ which diverges logarithmically with L 
as it was discussed above. These constants are actually essential for evaluating the 
PDF of the 2D-XY model [3] as they are related to the moments of the distribution. 
For example, magnetization goes slowly to zero with the system size like < m >= 
(NC)~ T ^ 87T ^ = e~ Tgi / 2 . The important point is that these constants depend only on 
dimension and lattice geometry. 

The scaling exponents of the spherical model in the range 2 < d < 4 are given by 
a = -(4 - d)/(d - 2), = 1/2, 7 = 2/(d -2),5 = {d + 2)/{d - 2), and v = l/(d - 2). 
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They satisfy the scaling relations a + 2/3 + 7 = 2, a = 2 — du, and 7 = v{2 — 77). In the 
critical region, we define therefore the rescaled parameters 



K — K c =: k^L- 1 ^ = k L-( d - 2 \ h =: hoL~^ v = h L^ d+2 ^ 2 , (12) 

where ko (which can be negative, in the paramagnetic region) and ho are external 
parameters. Similarly, h follows the same scaling as h, or h =: JiqL"^ 2 ^ 2 . Therefore 
the magnetization scales like < m >oc (K — oc h 1 ^ 6 oc L~( d ~ 2 )/ 2 . These arguments 
imply that the different moments and cumulants of order p scale like < m p >oc L~K d - 2 )/ 2 
in the critical region. Cumulants k p can be computed using the successive derivatives 
of $: 



1 ^ _. = T- P (d-2)/2 

NP-*dhr~ p 



(13) 



and the asymptotic distribution defined in the introduction can be identified with 



the function Q using a Fourier transform (see Appendix A) 



Q(9) :-. 



dX 
2^ 



exp 



Z—/ ~v 



{-i\y 



p>2 K 2 



p/2 



p\ 



(14) 



which is independent of the system size. We first evaluate the saddle point solution 
of Equation (JSJ), using the scaling Equation ( TT2|) and assuming z small: 



—T 



1 



2N ^0 d ~ Ea COS Qa 



+ 



q^O 



1 



2N^(d-Z a ^sq a ) 2 



+ 



1 



2Nz 



hi 



r-(d+2) + # + fc r-(d-2) + 



0. 



4iiT c z 2 " ■ — — - ( 15 ) 

A first analysis of the previous equation leads to the identification of the scaling 
law z ~ L~ 2 since for example z~ 2 L~^ dJr2 > ~ L-( d - 2 ) by comparison of terms on the 
second line. All other terms scale with the same exponent 2 — d. This means that 
d&/dz scales like L^^ 2 ^ and $ ~ L~ d at the critical point. Therefore the solution 
z(h) can be expanded as a series in the inverse power of L: z(h) = L^ 2 zq + • • • (and 
z(h) = L~ 2 z + ■ ■ ■ as well). We obtain at the leading order 0(L~( d ~ 2 ^) 



\gx + \z ~g 2 L-^ - \zl~g,L-^ + \zl~g,L~^ - • 
h 2 



1 

2~z~ 



L -(d-2) _ ^-(d-2) + ^ + fc L -(d-2) = , 

AK c zi 



(16) 



The zero th order term K c = g±/2 is canceled, and we obtain an implicit equation for 
z o — z o(ko, h ) and z = z (ko, h ) as function of the reduced temperature and magnetic 
field: 
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1 1 ~ 1 2~ 1 3~ h 2 

~2F + 2 Zm ~ 2^ 3 + 2^ 4 - ■ ■ • = ~ ^ (17) 

and a similar expansion for zq as well. We can express the previous equation in a 
more compact form, using the definition of the g p s and summing the series over p 

2z + 2 2^ i + 2oGqL - 2 AKc zl ^ ^ 

For the conjugate field of the absolute value \m\, we obtain the same expression 
with an additional term associated to the error function 

1 i z GlL-* _ hi hp ej^j) 

+ 4^ ° + 2^f 3 / 2 l + erf(^)' 1 ' 

When /i = 0, or h = 0, and k large, the previous equations give the asymptotic 
solutions zo ~ zo ~ l/(2/co). Numerically we solved Equation (fT8l) and Equation (fT9l) 
recursively to obtain the real solution z$ and zq. 

4. Cumulant expansion and asymptotic distribution 

In this section we derive the cumulants from Equation (|T3|) and given function $, using 
the constraint d^/dz = of the saddle point imposed for any conjugate field h or h. 
Indeed the identity 

(P(dM<h) t h,0)) = 

dhP 1 ' 

is valid for any integer p > 0, which simply means that at the saddle point value 
z = z(h), d z Q is always zero as function of h. In particular, we can show recursively in 



Appendix B that for p > 2 

1 d p $ 1 „ t <9 p $ 



K 



p A^" 1 d/i? A^p- 1 ^ p dz k dhP 

Knowing the function $, it is straightforward to obtain the cumulant expression by 
derivation and summation. It is worth noting that cumulants depend only on z{h) and 
its first derivative. From the scaling form of the distribution function, the scaling law 
for the cumulants (or equivalently for the moments) is k p ~ i j -p{ (1 - 2 )/ 2 . For example, 
Kx =< m >= h/(2Kz) = jk 1 ^^^ 2 ^ 2 = KiL~ (d ~ 2)/2 . We can generally use the fact 
that dz(h)/dh = dz (h )/dh x L^ 2 ^ 2 =: z' Q (h Q ) x L^ -2 )/ 2 , to calculate the regular part 
k p of the cumulants near the critical point, and obtain for example, using Equation (IB.5j) . 
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lp-2 



-Km 



pi 



2pz p 



+ 



4K c z 



p-i 



+ il v - 

2P ^ (z + L*G-iy 



from which we deduce the second cumulant 



u 'P— i 
2K c z p 



h 2 z' p 



(22) 



v'2 



^2 



2^ 



+ 



1 



/2 



2AT c z 

for order parameter m, and 



hoz' ^ K\z'q 



K c zl AK c z 3 



(23) 



?'2 



1 



4 



£ - | 1 | ° V — - 

2 2^ 2 2AT c z 2 ^ (zo + L 2 ^- 1 ) 2 tt^c^o V^ 3/ 



-3/2 



(24) 



for order parameter \m\, when h = for simplification. In the former case when 
ho = 0, we can show that z' (0) = 0. Indeed, deriving the saddle-point solution 
Equation ( |T8l) . we obtain 



h 2 

K c zl 



2 r-4 



+ 



z n 2 ^ ^ (1 + z G q L- 2 -> 



K c z 2 



(25) 



for conjugate field ho and 



-/ 

2n 



/^2 T -4 

^(l + *G,L-»)» 



3/2 ' 



(26) 



by deriving Equation ([19]) and taking /iq = afterward. The former equation 
gives z' (0) = 0, then only k.2 = 1/(2K c zq) is non-zero in this limit which implies that 
distribution Q(9) is purely Gaussian with a variance depending on an implicit equation 
for zo- However in the latter case Equation (1261) . z' Q (0) ^ instead and Q{6) is non- 
Gaussian since all the moments do not vanish. We can in general derive formally Q(6) 
from definition Equation (1A.3|) . using the expression Equation (121]) for k p , and the 
characteristic function. We also redefine the free energy <3>o := L d &(zoL~ 2 , hoL~( d+2 ^ 2 ), 
with 



^o(z ,h ) = -^log2;o + ^-^2 



H - log 1 + 



L 2 



L 2 



+ 



h 2 



AK c z 



+ k z + est. (27) 



For the distribution of \m\ and conjugate field ho, we find instead 
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$0(^0,^0) 



1 , 1 

-2 1Og " + 2 



+ log 



1 + erf 



h 



— log 1 + 



V- 



L 2 



+ 



hi 



AK c z 



+ est. 



(28) 



The characteristic function H appearing in the definition of Q(9) in Equation flA.lj) 
has the scaling form H depending on $ ( or ^0) an d its derivatives 



p>2 K 2 P- p>2 K 2 P- k=0 



(29) 



The series over p can be performed directly if we use for example the Fourier 
transform of $0 



$0(^0,^0) 





dqdtu ~ . . .. . , , 
— — ~®o{q, u) exp(iqz + iuh ), 
(2vr) 2 

dzodh $ (z , h Q ) exp(—iqz — iuh ). 



(30) 



Indeed, we obtain formally the functional expression of the characteristic function 
depending on $ only 



2^ -P/2 , k p 

p>2 K 2 P l 




dqduj ~ 



j^Mq, w) exp(iqz + iuh ) £ T^T" C k p {iq) k {iuf~ k z' Q {h) k 
dqdu ~ 

W) " 




p>2 K 2 P- k=0 

-i\y 



$o(q,uj)exp(iqz + iujho)^2—^ — (iqz' + 



p/2 , 

p>2 «2 P l 



(31) 



<I>o( -u - z-^S=, h - z-^=) - $0(^0, ho) + i^= [z'o^- + 7^- ) *n( -n- /'n ) 

V K 2 V K 2 



'«2 



This relation satisfies the constraint i?(0) = coming from the normalization of 
the PDF. It is also straightforward from Equation (1A.3|) to obtain the general form of 
the limiting distribution for conjugate fields ho and ho as well by replacing $ by $0 



Q{0)= J ^exp 



iX<0 
Xz' 







+ $0(2:0 -i^=,h -% 



% dh 
X 4 



$0(^0,^0 
- $0(^0,^0) 



(32) 



'K 2 V K 2 

This result is convenient since the Fourier transform depends functionally only on 
$0 which can be written for different lattice geometries or coupling distributions between 
the nearest-neighbor spins, by only modifying the structure of the Green function G q 
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Figure 2. (Color online) plot of Zq, z' , k 2 , and asymptotic slope u* — y/K^z /z' Q as 
function of fco for the distribution of |m|. In the limit of large fcoi u* approaches the 
constant u*^ — \/ir(l — 3/ir)/2 ~ 0.266 (dashed orange line, see text). 



for example. Also only zq and z' (or zq and z' ) are necessary for the evaluation of 
this expression by solving Equation (jl~8j) and Equation (|25|) (or Equation (jl~9j) and 
Equation (|26|) ). The distribution is centered around the value 



Co 



d 



%dh 



(33) 



which is equal to c = h /(2y/K^K c z ) using Equation ( 1271) or c = 1/ ' ^k^K^zq 
using Equation (1281) when h = for simplification. Typical behavior of quantities z Q , 



and «2 are plotted in Figure |2] after solving Equation ( fl~9l) recursively. 



5. Asymptotic limits and numerical results 

The asymptotic analysis of Equation ( )32|) when #;^>lor#^— lis useful to determinate 
the behavior of the PDF for large deviations and obtain its universal characteristics. 
The standard method is to derive the argument in the exponential of the integral 
Equation fl32l) with respect with A and look for the dominant saddle point in the complex 
plane. This is equivalent to analyze the existence of solutions of the following equation: 



6 + 



1 d_ 

/n^dh^ 



$o(zo,^o) 



+ 



/k2 dz 
1 d% 
fW 2 dh 



V*0 



hr 



hn - i- 



X 



/«2 



(34) 
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The detailed analysis of this equation is given in Appendix C We find two kinds 
of behavior. For 9-^—1 the PDF is exponentially decreasing Q(9) ~ exp(u*9) with 
coefficient u* (and u*) equal to u* := v^W^o > (u* := \^z /z > 0)- The saddle 
point structure can be checked numerically by fitting the asymptotic function given 
below by Equation (I35p where corrections have to be taken into account for moderate 
9 < in addition to the dominant exponential term. In particular three dominant 
contributions are present corresponding respectively to the linear behavior discussed 
above, plus a term in ^\9 + cq\, and a logarithm log \ 9 + cq|: 

Q{9) ~ exp (u*9 + ~jL b^^\9 + c | - \ log \9 + c \) , 9 < -1. (35) 

In the opposite limit, when 9 ^> 1, the PDF falls off exponentially with a stretched 
exponent equal to d/(d — 2) = 3 in three dimensions. The saddle point evaluation in 



this case is detailed in Appendix C 



, pz- v d/(d-2) 

Q{9) ~ exp (- /3(d) (^^J ),#»!, (36) 

with eventual corrections. Here the expression of coefficient /3(d) is given by 
Equation flC.111) . and is simply equal to 0(3) = 8n 2 in three dimensions. These results 
can be checked numerically after plotting the Fourier integral Equation (|32|) . 

In Figure [3] is represented the distribution for m and for two different sets of 
parameters (ho,ko). Q(9) is typically non-Gaussian, with a exponential behavior for 
negative deviations, as expected, and more pronounced falloff for positive deviations 
for which a stretched exponential with a cubic exponent was performed adequately. 
The exponential behavior is typical to extreme value statistics when one studies the 
extrema distribution of a set of independent random variables, such as the Gumbel 
distribution. This distribution has a more pronounced asymptotic double-exponential 
falloff with a coefficient which is integer. Here the coefficient of 9 3 is non-integer and 
depends on the saddle point structure. When field ho is increased, negative deviations 
are enhanced, whereas the curve falls off more rapidly for 9 > as the coefficient 
(v^/^o) 3 m Equation ( |36l) increases. 

The distribution for \m\ follows the same trend. For large and negative deviations, 
the saddle point solution is given by Equation ( 1C.6I) . with an asymptotic behavior also 



corresponding to Equation (1351) and Equation (136|) . Figure H] represents typical examples 
for k = 0.1 and ko — 1 at zero field h = 0. The PDF is clearly non-Gaussian, and, 
for deviations left to the most probable value, two distinct behaviors separated by a 
crossover interval. When 9 is largely negative, the slope of the exponential asymptote is 
given by u*, whereas there exists a plateau-like regime where the distribution, at least 
for k > 1/2, takes significant values in a whole range of 9 values. However, when fc 
decreases, the asymptotic regime is reached only for negative deviations that become 
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10 n 



~i — i — i — i — i — i — i — i — i — i — i — i — i — i — i — i — i — r 



i — i — i — i — r 





h =l,k =0.1 
h =0.1,k n =0.1 



-2 







Figure 3. (Color online) plot of distribution Q(9) Equation (|32|) for m and for 
two different sets of parameters (ho,ko) (L — 10). Dashed lines correspond to 
approximation curves given by asymptotic regime detailed in |Appcndix C For negative 
deviations, the behavior is found to be exponential Equation Q35[) with u* ~ 7.396 for 
(ho = 0.1, fco = 0.1) and u* ~ 1.44 for (ho = 1, fco = 0.1), close to expected values 
u* = 8.15 and u* = 1.45 respectively. For positive deviations 9 > 1, curves are fitted 
with a stretched exponential with cubic exponent Equation (|36[) . 



larger and larger, as it is shown in the set of curves Figure [5] where we have varied fco 
almost continuously. 

In particular, in this figure, the limit A; ^> 1 is seen to be reached for k values 
close to 2, and the limiting distribution takes a universal form. 

Indeed the saddle point values when h Q = can be computed exactly in this limit 
since z ~ l/(2fc ), z' ~ a/Io, and k 2 ~ 1/%, and u* approaches the numerical constant 
u*^ = y/ir(l - 3/tt)/2 ~ 0.266, see |Appendix C\ Equation (ICTfj) . We find that the 
limiting value Qoo of the distribution is given by 



Qoo{9) 



2tt 
ttA 2 
4(1 -iX) 



exp 



-l=)-llog(l-,A 



+ log 



erf 



i^/nX 



(37) 



• 2^1 -i\< 

This integral, independent of ko, is a universal function of 9 with numerical 
factors only. It is plotted in Figure |6] and presents a dominant contribution in the 
interval —2.1 < 9 < 0.9, with a sharp decreasing behavior outside this interval. The 
plot suggests that there may exist a cutt-off at 9 = 9 C ~ 0.86 above which the 
distribution vanishes. We have approximated quite accurately this function by the 
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Figure 4. (Color online) plot of distribution Q(9) Equation (I3"2"j) for |m| and for 
two different parameters fco (ho = and L = 10). Dashed lines correspond to 
approximation curves given by the asymptotic regime. For negative deviations, the 
fits are given by the exponential behavior Equation (1551) with u* ~ 0.26 for k = 1, 
close to the expected value u* — 0.279. For fc = 0.1, the asymptotic regime is not 
reached since we expect a slope u* = 0.62, lower than the slope found in the interval 
[—3, — 1]. For positive deviations 6 > 1, curves are fitted accurately with a stretched 
exponential with cubic exponent Equation (|36[) . 



Ansatz Qoo(0) w exp{a6> - 6/(0 - C )}/(9 C - 6) a in the interval 9 > -2, with fitted 
parameters a = 0.055, b = 0.0351, # c = 0.8533 and a = 0.826. 

Finally, in Figure [TJ we have plotted the distribution of m for a series of several ho 
values at k Q = 0. As discussed above, the distribution for (ho = 0, k Q = 0) is a Gaussian, 
and differs from the Gaussian form as the field is increased, and presents asymptotically 
the exponential form Equation (I35j) for large negative deviations. However the structure 
is different from Figure [5] since no plateau regions are present. 

6. Conclusion 

In this paper, we presented an analytical method to compute the distribution function 
of the magnetization in the spherical model at criticality, with a numerical application 
in three dimensions. The advantage of this spin-correlated model is that the critical 
region can be explored analytically leading to important properties of the PDF. In 
particular, the exponential behavior for large and negative deviations are characteristic 
of the presence of soft modes |q| <C 1 destroying the long-range order as seen in the 2D- 
XY model, with additional corrections not seen usually in extreme statistics distributions 
of uncorrelated variables. For larger deviations, the PDF is falling off more rapidly than 
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Figure 5. (Color online) multiple plots of the distribution Q(9) for \m\ and for several 
values of fco (ho = and L = 10). A crossover occurs around k — 0.5 below which the 
plateau-like feature of the distribution is smoothed out. The large ko limit is almost 
reached when fco is larger than 2 (see also Figure [5]) . 




Figure 6. (Color online) plot of the distribution Q oc (9) Equation (1571) in the limit 
of large fco (ho = and L = 10). Dashed blue line corresponds to the approximate 
function Q 00 (6') s» exp{a6>— b/ (0—9 c )}/(9 c — 6) a in the interval 9 < 9 C , with parameters 
a = 0.055, b = 0.0351, 9 C = 0.8533, and a = 0.826. 
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Figure 7. (Color online) multiple plots of the m distribution Q(9) for several values 
of ho ( fco = and L = 10). The distribution at ho = is purely Gaussian 
Q(9) = exp(~6» 2 /2)/\/2^ (see text). 



a Gaussian, as the long-range order is more difficult to restore. The stretched exponent 
of the exponential falloff depends on the dimension like d/(d — 2). When d = 4 the 
Gaussian behavior is recovered, and when d is close to two dimensions, the exponent is 
diverging, and we expect formally a double exponential falloff in this limit for positive 
deviations, as for the 2D-XY model and 2D Ising model as well [23]. By inspection of the 
partition function in the Fourier space Equation ([5]), we can reformulate the problem and 
associate the PDF as the distribution of a set of uncorrelated Fourier modes S" q which 
appear in the integrand. However these modes are coupled by the variable of integration 
s coming from the constraint over the spins, which leads, a priori, to a PDF which is not 
a standard distribution. It seems nonetheless that using a saddle point analysis in the 
large size limit we obtain the PDF of quasi-independent modes not evenly distributed, 
which implies strongly non-Gaussian characteristics in the general case [19J. We found 
for the spherical model that the PDF belongs to a class of functions closely related to 
Gumbel but with a stretched exponential for positive deviations instead of a double 
exponential form when the dimension is larger than two. Moreover the corrective terms 
found asymptotically for negative deviations in Equation ( |35l) may originate from the 
correlations between modes and induced by the constraint over the spins. 

The main result Equation ( 1321) . giving formally the Fourier expression of the PDF, 
can be applied for any model where a saddle point analysis is exact, since the cumulants 
can be expressed simply with a generating function, see Appendix B This model can 
moreover be modified in order to incorporate additional couplings or crystal fields for 
example. In this case, only the rescaled free energies given by Equation ( |27|) and 
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Equation (1281) have to be reformulated to accommodate the modifications made in the 
original Hamiltonian. 

Appendix A. Limiting distribution 

Cumulant /%, are defined by the characteristic function H of instantaneous spin order 
parameter ^^Si/N 

E [e iX ^ N ] = exp[#(A)] = exp ( $> p ^ ) (A.l) 

\p>i p ' J 

where E[.] is the thermal average operator over the different spin configurations 
and k p the cumulants. The distribution P(m, L) can then be expressed as the Fourier 
transform of exp [if (—A)] 



P(m,L) :-- 



exp[zAm + H(— A)] 



dX 
2^ 

dX 

^ exp 



i Am 



p>i 



-ixy 



pi 



(A.2) 



If k p = for all p > 2, P(m, L) is Gaussian and K2 is the variance a 2 



(m—Ki) 



J /y/fy, where Q{9) 



P(m, L) = exp \ - {m ~^' j /V%™2 or P(m, L) = Q [6 := 

exp(— 9 2 /2) j \f}/K is the normalized Gaussian distribution function. In general, we expect 
the distribution to be non-Gaussian if k p > 3 ^ 0, and P(m, L) can be put generally into 
the following form using the rescaled cumulants defined in Equation (Fl3l) 



P(m,L) = 



dX 



dX 



== exp 



iX 



[m — Ki) 



+ £5 



-ixy 



pi* p! 

p>2 ^2 ^ 



7= ex P 



K P (-ixy 



p>2 K 2 F 



with /3/v = (d— 2)/2 for the spherical model. 



V^2 



(A.3) 



Appendix B. Cumulant expression 

In this appendix, we derive the exact differential relation Equation (1211) . Consider thep th 
cumulant of the distribution as the derivative of the (p — l) th one, using the assumption 
that Equation (12ip is correct for k p _i 
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\fe=o / 

, i k ^(h) k+1 + a ,rS , ^+ ,1 , r kz'(h) k - 1 z"(h) . 

dz k+1 dhP- 1 - k v ' dz k dhP~ k v y dz k dhP~ 1 ~ k v 7 



fc=0 



The first two terms in the last line can be brought together using the binomial 
formula Czi\ + Cp_i = Cp arid rearranging the summation index k. We then obtain 

p f)P(b p ~ 1 c)P~ l & 

N p -\ = E + A*) E ^i^^^W- 1 - (B.2) 

fc=0 fc=0 

Using kCp_ t — (p — l)CpZ2i we can show that the last term is equal to 

fe=i fc=i 

= A/0(p-i)^3t- (B.3) 

This term is equal to zero at the saddle point value, proving the recurrence for 
Equation ( 12TT) at order p. Now, if we consider the first expression for $ in Equation ((7j), 
the dependence in h is quadratic, and only 3 terms remain from Equation f|2T|) 



Kp 



NP~ 



/p gPg /p _! 1 gPg /p _ 2 



• (B.4) 



Using for example the exact expression for $(2, /i, 0), we can compute the cumulants 
as function of the saddle point z and its derivative z'ih) only 

(-1)p z'p z'p~ 2 z'p 1 ^ 1 

p\ ^ p ~ 2 P Np Z p + AKNp- 1 zp- 1 + 2p~NP ^ (z + d - J2 a cosq a ) p 

hz'P- 1 h 2 z' p , . 

+ 71^=T-^T - B.5 



2KNp~ 1 zp AKNp~ 1 zp+ 1 ' 
This result can be used in Equation ( 1A.2|) to obtain the characteristic function 



H{— A) by re-summation over p of the cumulant series. 
Appendix C. Asymptotic analysis 

The asymptotic solutions for the saddle point equation Equation ( 1341) can be derived 
when 9 <C — 1 and 6 ^> 1. In both cases, the left hand side of Equation fl34|) proportional 
to # is diverging negatively and positively respectively, which leads us to look for a 
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diverging solution for the right hand side as well. For the two cases considered in this 
paper, Equation ( 1511) is equal to 



and 



6 + 



1 



>K2 2K C Z 



+ 



4 


1 




2Z, 


1 


H 



ik 2 2K C Z ' 



1 



Hi 



2 U L4 1 + Z ° GJL2 



+ k 



(CI) 



+ 



y/K 2 irK c Z 

Ho 



'«2 



2Z. 



-H$/(4K C Z ) 



2^K c Zl' 2 1 + erf 
where 



Ho 



2y/K c Z 



Gl 



Zo 



2 ^ L 4 1 + Z G q /L 2 4AT C Z, 



o 



+ 



'«2 



H n 



2K C Z 



1 



-i? 2 /(4X c Z ) 



WcZol + erf — ^ 



(C.2) 



2y/K c Z J _ 







. A 






= z 


Zo 


Ho 


= h 




Ho 



zo 



'«2 



A 



K>2 



(C.3) 



Considering first the regime 9 <C — 1, let assume in the following that zo and are 
positive (as well as z Q and z' ), which will be realized numerically by solving Equation f|T8|) 
and Equation (12"oT) . We are looking at diverging terms in the previous saddle point 
equation coming from the inverse powers of Zq — or Zo — for example. Considering 
the integration path C~ in the complex plane for A (see Figure ICTI) . it is clear there 
is a special point on the negative imaginary axis A = — iu* with u* := v^W^o > 
(w* := ^/k^Zo/z'q > 0) for which Z = and the right hand side of Equation ( 1C.1I) (and 
Equation (1C.2j) ) is singular. 

The most singular term appears to be the one proportional to 1/Zq oc l/(w* — u) 2 , 
after setting A =: —iu with u close enough to u* , and which gives the estimate for the 
saddle point solution 



u 



u 



u*\ho — M*/v^2| 



(C.4) 



2y/K c z \e + co\ 

Replacing this value in Equation ( )32|) and integrating the Gaussian fluctuations 
around the saddle point, we obtain the dominant contribution for the large negative 
deviations Equation ( 1351) . 

The asymptotic behavior of the PDF is therefore exponentially decreasing, with a 
coefficient equal to u* . The same analysis for Equation (1C.2I) leads to a similar result, 
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Figure CI. (Color online) modified path of integration for the saddle point analysis. 
C + is the path chosen when 9 >• 1 and when 9 <C — 1. 



with coefficient u* instead. The inverse powers in Equation (IC.fj) and Equation (IC.2j) 
appear to be indeed only in l/Z , 1/Z^ 2 , 1/Z^ 2 , and 1/Zq, with the latter one being 
the dominant contribution. The additional two terms in Equation (1C.2j) proportional to 



-H$/(4K C Z ) 



-u* 3 /[4K c k 2 z (u* -«)] 



1 + erf 



H 



1 - erf (u* 3 / 2 /[2^K c k 2 z {u* - u)} 



(C.5) 



,2y/K c Z , 

are singular in the limit since both numerator and denominator vanish. 

However this ratio can be evaluated using the asymptotic behavior of the error function 

— 1 /2 

for large and negative argument. This ratio diverges like 1/Z and elevates the power 



of the factor 1/Z^ 2 in the second line in Equation (1C.2|) . giving rise to a dominant 



contribution proportional to 1/Zq, in addition to the identical contribution proportional 
to Hq/Zq found in the first line. Finally, the saddle point solution is given by 



u ~ u* — 



2K c zo\9 + Co\ 



(C.6) 



instead of Equation (10.4[) and for h = only. Since Z is small, the sum term over 
modes q in the saddle point equations is regular since it behaves in this limit like Z g 2 
with g 2 finite. It is no more the case in d > 4 where g 2 is diverging logarithmically with 
the system size and where a careful different analysis has to be made for the new saddle 
point, and where one expects to find a Gaussian behavior with a 8 2 contribution. 

Considering the PDF of \m\ in absence of field h = 0, we can estimate the 
asymptotic value of u* when ko is large. Indeed, from Equation (|T9l) we obtain the 
approximation zq ~ 1/(2A; ) 1, and from Equation ( 1261) . z' ~ a/ z /(tcK c ). In 
addition, the second cumulant is equal to k 2 
the universal limit 



'1 — 3/7r)/(2K c z ). Then we obtain 
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lim u* =: = J - ( 1 - - ) ~ 0.266. (C.7) 



fco->oo v V 2 \ 7T/ 

In the regime of large and positive deviations 9 1, the previous analysis can 
not be valid anymore because of the sign of 9. Instead we are looking for a path 
C + in the upper complex plane 3(A) > where Z is real and large. Indeed since 
9 is large, we can assume that the modulus of A =: iu becomes large as well, so 
that Zq ~ uZq/t/kq ^> 1. In this case, terms proportional to inverse powers of Z$ 
in Equation (IC.ip and Equation 0U.2p are small or finite, except the sum over the q 
modes. This sum contributes like 2~^ q ^o ^q/ ^ 2 which diverges with the system size when 
d > 2. Indeed g± = J^q^o^q/-^ ^ s finite for d > 2 and diverges logarithmically when 
d = 2. Therefore the former sum is divergent with Z . Since Z > when A = iu with 
u > 0, the sign of the function 

^H£§T7z|^ (a8) 

is positive and consistent with the sign of 9 on the left hand side of Equation (IC.lj) , 
and the saddle point equations reduces to 9 ~ z' Q F(Z ) / '(2\/k2). The divergent part 
of F(Z ) can be evaluated exactly in the continuous limit. Indeed, we can rewrite 
Equation (1C.8|) as an integral 



nz a ) - ^ f * (c.9) 



WJvt/L q i LH + 2Z /( q ?L*y 

where = 2n d ^ 2 /T(d/2) is the hyper-spherical volume. Performing the change of 
variable qL/^2Z — > g and taking afterward the limit L — > oo, we obtain typically 



^ - wL^^ 1 " - PiW^CT - (C10) 



where the integral can be computed in general dimension 2 < d < 4 after 
substituting the lower bound with zero. Then the saddle point is given by 



{(2/7r)(27r) d / 2 r(d/2)cos[7r(d + l)/2]} i^j 9 2 '^ » 1 

0(d) (V^) a/(*-2). (C.ll) 



After inserting this value in Equation f )32|) . we finally obtain the dominant behavior 
of the PDF given by the stretched exponential Equation (I36p . 
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